/*
 * File:   fe_matr.h
 * Author: boyan
 *
 * Created on December 11, 2007, 12:29 PM
 */

#ifndef _FE_MATR_H
#define	_FE_MATR_H

#include "glodef.h"
#include "grid3d.h"
#include "vecs.h"
#include "fe.h"

/* different entry operations during matrix assembly */
enum {
    eo_add, /* the standard */
    eo_set, /* for the df-basis we need the set sometimes */
    eo_ignore,
    eo_bc, /* set diagonal el to 1 and off-diag in the same row and col to 0 */
    eo_condense
};

typedef void get_eo_t(const grid3d_t *grid, unsigned eltype, const int *nodes, unsigned **eos);
get_eo_t eo_bc_dir, eo_add_all;
get_eo_t eo_condense_bc_dir, eo_condense_add_all;

typedef struct _fe_matr_t * fe_matr_p;
unsigned get_el_type(fe_matr_p M);

#ifdef	__cplusplus
extern "C" {
#endif

    fe_matr_p fe_matr_alloc(unsigned type, local_matr_t *lm, get_eo_t *get_eo,
            grid3d_t *grid, int symmetric, char *mattyp);
    void fe_matr_free(fe_matr_p matr);
    //    void fe_matr_structure(const grid3d_t *grid, int nmatr, fe_matr_p *matr);
    void fe_matr_assemble(const grid3d_t *grid, int nmatr, fe_matr_p *matr);

    int fe_matr_solve_vector(const fe_matr_p matr, vector X, vector F, vector T, double TOL);
    int fe_matr_solve_scalar(const fe_matr_p matr, scalar X, scalar F, scalar T, double TOL);
    /* Solves matr*X=F.  The boundary values of X are assumed to be in F.
     * X is used as initial guess without any modification.
     * T is auxilliary, TOL is passed to sparse_cgm or sparse_bcgm.
     */

    void fe_matr_print(const fe_matr_p matr, const char *fname, const char *fname_bc);

    void project_grad_vec(const grid3d_t *grid, unsigned type, const fe_matr_p mass,
            vector u, vector gu0, vector gu1, vector gu2, vector w, vector t);
    void project_grad_scal(const grid3d_t *grid, unsigned type, const fe_matr_p mass,
            scalar p, vector gp, vector w, vector t);

#ifdef	__cplusplus
}
#endif

#endif	/* _FE_MATR_H */

